library(topGO)
library(stats)
library(knitr)


setwd('/Users/mmayers/projects/n15_mice/data/')

data.dir <- '../data'

gene2GOID <- readMappings('clusterID2GO.map')
unenr_pvals <- read.csv(file.path(data.dir, 'unenriched_pvals.csv'))
unenr_pvals$log2FoldChange <- log2(unenr_pvals$ratio)
unenr_pvals$padj <- p.adjust(unenr_pvals$p_value, method="BH")
GO2geneMF <- annFUN.gene2GO('MF', gene2GO=gene2GOID)
GO2geneBP <- annFUN.gene2GO('BP', gene2GO=gene2GOID)

do some stuff with fold changes

unenr_pvals$adjFC <- unenr_pvals$log2FoldChange
n <- sum(unenr_pvals$padj > .05 & abs(unenr_pvals$log2FoldChange) > 2)
unenr_pvals[unenr_pvals$padj > .05 & abs(unenr_pvals$log2FoldChange) > 2, ]$adjFC <- rep(0, n)
#set.seed(2222)
#unenr_pvals[unenr_pvals$padj > .05 & abs(unenr_pvals$log2FoldChange) > 2, ]$adjFC <- runif(n, -1, 1)
hist(unenr_pvals$log2FoldChange, 15)

hist(unenr_pvals$adjFC, 15)

Metadata on Cluster Info - Use this for filtering of certain types of proteins

sampleTable <- read.csv(file.path(data.dir, 'filt_metadata.csv'), row.names=1)
sampleTable$technical <- as.logical(sampleTable$technical)
locusTable <- read.csv(file.path(data.dir, 'loci_annot.csv'))

Look at Bacterial Proteins only

Mouse Human subset seems to be a bust, can’t figure out good way to account for tiny background)

tmp <- locusTable$mouse_human == 'False' & locusTable$lca != '35823' # Filtering out anthrospira differences 35823
tmp[is.na(tmp)] <- TRUE
bac.prots <- locusTable[tmp, ]$X
bac.pvals <- unenr_pvals[unenr_pvals$X %in% bac.prots, ]
geneList.bac <- bac.pvals$padj
geneList.bac <- bac.pvals$adjFC
names(geneList.bac) <- bac.pvals$X

topDiffGenes <- function(score) {return(score >= 2)} # scoring function, anything with > 2 logfc is sig anyway

TCell vs RAG - Unenriched

geneList.bac <- -1*geneList.bac

MF first

ontology <- 'MF'
GOdata.bac <- new('topGOdata', ontology = ontology, allGenes = geneList.bac, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.ks <- runTest(GOdata.bac, statistic = 'ks')
a <- GenTable(GOdata.bac, result.bac.ks, topNodes = 10)

kable(a)

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0004619	phosphoglycerate mutase activity	42	0	4.77	1.4e-26
GO:0000287	magnesium ion binding	46	0	5.23	8.6e-09
GO:0005506	iron ion binding	80	0	9.09	7.2e-07
GO:0016491	oxidoreductase activity	291	15	33.07	1.3e-06
GO:0008184	glycogen phosphorylase activity	14	0	1.59	1.6e-06
GO:0004634	phosphopyruvate hydratase activity	34	0	3.86	5.0e-06
GO:0046933	proton-transporting ATP synthase activit…	14	1	1.59	4.1e-05
GO:0016639	oxidoreductase activity, acting on the C…	29	0	3.30	4.3e-05
GO:0004618	phosphoglycerate kinase activity	27	3	3.07	6.3e-05
GO:0003677	DNA binding	61	0	6.93	6.8e-05

for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac, term, ranks=TRUE))
}

## [1] "phosphoglycerate mutase activity"

## [1] "magnesium ion binding"

## [1] "iron ion binding"

## [1] "oxidoreductase activity"

## [1] "glycogen phosphorylase activity"

## [1] "phosphopyruvate hydratase activity"

## [1] "proton-transporting ATP synthase activit..."

## [1] "oxidoreductase activity, acting on the C..."

## [1] "phosphoglycerate kinase activity"

## [1] "DNA binding"

#showGroupDensity(GOdata.bac, "GO:0008184", ranks=TRUE) #Glycogen phosphorlase
#showGroupDensity(GOdata.bac, "GO:0004618", ranks=TRUE) #phosphoglycerate kinase

ontology <- 'BP'

BP - Tcell vs RAG: Unenriched

GOdata.bac <- new('topGOdata', ontology = ontology, allGenes = geneList.bac, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.ks <- runTest(GOdata.bac, statistic = 'ks')
a <- GenTable(GOdata.bac, result.bac.ks, topNodes = 10)

kable(a)

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0006096	glycolytic process	148	25	15.31	2.3e-18
GO:0005975	carbohydrate metabolic process	320	48	33.11	4.0e-14
GO:0008152	metabolic process	1017	92	105.24	2.1e-06
GO:0055114	oxidation-reduction process	312	28	32.29	1.8e-05
GO:0015986	ATP synthesis coupled proton transport	15	1	1.55	0.00028
GO:0006520	cellular amino acid metabolic process	78	0	8.07	0.00031
GO:0006108	malate metabolic process	6	0	0.62	0.00037
GO:0008654	phospholipid biosynthetic process	13	0	1.35	0.00072
GO:0006021	inositol biosynthetic process	13	0	1.35	0.00072
GO:0006351	transcription, DNA-templated	63	0	6.52	0.00077

for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac, term, ranks=TRUE))
}

## [1] "glycolytic process"

## [1] "carbohydrate metabolic process"

## [1] "metabolic process"

## [1] "oxidation-reduction process"

## [1] "ATP synthesis coupled proton transport"

## [1] "cellular amino acid metabolic process"

## [1] "malate metabolic process"

## [1] "phospholipid biosynthetic process"

## [1] "inositol biosynthetic process"

## [1] "transcription, DNA-templated"

#showGroupDensity(GOdata.bac, "GO:0006096", ranks=TRUE) # Glycolitc Process
#showGroupDensity(GOdata.bac, "GO:0005975", ranks=TRUE) # Carbohydydrate metab
# Mostly Carb metabolism?  lipid biosynthesis?Nah, those look crappy

RAG vs Tcell - Unenriched

geneList.bac <- -1*geneList.bac

MF

ontology <- 'MF'
GOdata.bac <- new('topGOdata', ontology = ontology, allGenes = geneList.bac, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.ks <- runTest(GOdata.bac, statistic = 'ks')
a <- GenTable(GOdata.bac, result.bac.ks, topNodes = 10)

kable(a)

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0005524	ATP binding	233	2	4.38	1.7e-24
GO:0003872	6-phosphofructokinase activity	25	0	0.47	1.8e-14
GO:0004478	methionine adenosyltransferase activity	29	0	0.55	8.1e-13
GO:0005198	structural molecule activity	151	0	2.84	5.5e-12
GO:0050242	pyruvate, phosphate dikinase activity	74	0	1.39	6.9e-12
GO:0016887	ATPase activity	34	0	0.64	2.5e-11
GO:0008878	glucose-1-phosphate adenylyltransferase …	15	0	0.28	7.0e-11
GO:0016301	kinase activity	128	1	2.41	5.3e-09
GO:0003924	GTPase activity	77	1	1.45	4.1e-08
GO:0003824	catalytic activity	1023	22	19.24	1.7e-06

for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac, term, ranks=TRUE))
}

## [1] "ATP binding"

## [1] "6-phosphofructokinase activity"

## [1] "methionine adenosyltransferase activity"

## [1] "structural molecule activity"

## [1] "pyruvate, phosphate dikinase activity"

## [1] "ATPase activity"

## [1] "glucose-1-phosphate adenylyltransferase ..."

## [1] "kinase activity"

## [1] "GTPase activity"

## [1] "catalytic activity"

BP - RAG vs Tcell: Unenriched

ontology <- 'BP'
GOdata.bac <- new('topGOdata', ontology = ontology, allGenes = geneList.bac, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.ks <- runTest(GOdata.bac, statistic = 'ks')
a <- GenTable(GOdata.bac, result.bac.ks, topNodes = 10)

kable(a)

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0006556	S-adenosylmethionine biosynthetic proces…	29	0	0.48	6.8e-14
GO:0071973	bacterial-type flagellum-dependent cell …	95	0	1.57	3.8e-13
GO:0006002	fructose 6-phosphate metabolic process	22	0	0.36	7.4e-13
GO:0006090	pyruvate metabolic process	222	0	3.67	1.1e-12
GO:0016310	phosphorylation	229	1	3.78	1.2e-12
GO:0005978	glycogen biosynthetic process	15	0	0.25	4.1e-11
GO:0009401	phosphoenolpyruvate-dependent sugar phos…	13	0	0.21	4.6e-07
GO:0044267	cellular protein metabolic process	130	2	2.15	1.1e-06
GO:0006810	transport	61	0	1.01	4.0e-06
GO:0055114	oxidation-reduction process	312	4	5.15	4.6e-05

for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac, term, ranks=TRUE))
}

## [1] "S-adenosylmethionine biosynthetic proces..."

## [1] "bacterial-type flagellum-dependent cell ..."

## [1] "fructose 6-phosphate metabolic process"

## [1] "pyruvate metabolic process"

## [1] "phosphorylation"

## [1] "glycogen biosynthetic process"

## [1] "phosphoenolpyruvate-dependent sugar phos..."

## [1] "cellular protein metabolic process"

## [1] "transport"

## [1] "oxidation-reduction process"

BioGlyCMK Rag vs Tcell

enr_pvals <- read.csv(file.path(data.dir, 'enriched_pvals.csv'))
enr_pvals$log2FoldChange <- log2(enr_pvals$ratio)
enr_pvals$padj <- p.adjust(enr_pvals$p_value, method="BH")

Fix fold change of insignifcant proteins

enr_pvals$adjFC <- enr_pvals$log2FoldChange
n <- sum(enr_pvals$padj > .05 & abs(enr_pvals$log2FoldChange) > 2)
#set.seed(1234)
#enr_pvals[enr_pvals$padj > .05 & abs(enr_pvals$log2FoldChange) > 2, ]$adjFC <- runif(n, -1, 1)
enr_pvals[enr_pvals$padj > .05 & abs(enr_pvals$log2FoldChange) > 2, ]$adjFC <- rep(0, n)
hist(enr_pvals$log2FoldChange, 15)

hist(enr_pvals$adjFC, 15)

Bacterial proteins only

bac.enr.pvals <- enr_pvals[enr_pvals$X %in% bac.prots, ]
#geneList.bac <- bac.pvals$padj
geneList.bac.enr <- bac.enr.pvals$adjFC
names(geneList.bac.enr) <- bac.enr.pvals$X

Tcell vs RAG - Enriched

geneList.bac.enr <- -1*geneList.bac.enr # KS goes from - to +
                                        # do this to put upreg in Tcell first

MF

ontology <- 'MF'
GOdata.bac.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.bac.enr, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.enr.ks <- runTest(GOdata.bac.enr, statistic = 'ks')
a <- GenTable(GOdata.bac.enr, result.bac.enr.ks, topNodes = 10)

kable(a)

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0016639	oxidoreductase activity, acting on the C…	70	0	0.45	2.5e-07
GO:0016491	oxidoreductase activity	367	2	2.36	6.1e-07
GO:0016829	lyase activity	160	1	1.03	1.2e-06
GO:0051287	NAD binding	22	0	0.14	5.1e-06
GO:0016462	pyrophosphatase activity	50	0	0.32	6.6e-06
GO:0050661	NADP binding	20	0	0.13	9.4e-06
GO:0005515	protein binding	16	0	0.10	0.00020
GO:0005215	transporter activity	25	0	0.16	0.00025
GO:0004197	cysteine-type endopeptidase activity	7	0	0.04	0.00135
GO:0005506	iron ion binding	52	0	0.33	0.00155

for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac.enr, term, ranks=TRUE))
}

## [1] "oxidoreductase activity, acting on the C..."

## [1] "oxidoreductase activity"

## [1] "lyase activity"

## [1] "NAD binding"

## [1] "pyrophosphatase activity"

## [1] "NADP binding"

## [1] "protein binding"

## [1] "transporter activity"

## [1] "cysteine-type endopeptidase activity"

## [1] "iron ion binding"

BP Tcell vs RAG: BioGlyCMK

ontology <- 'BP'
GOdata.bac.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.bac.enr, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.enr.ks <- runTest(GOdata.bac.enr, statistic = 'ks')
a <- GenTable(GOdata.bac.enr, result.bac.enr.ks, topNodes = 10)

kable(a)

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0006520	cellular amino acid metabolic process	140	0	0.65	7.7e-13
GO:0055114	oxidation-reduction process	369	2	1.72	7.6e-10
GO:0006508	proteolysis	19	0	0.09	7.4e-07
GO:0006814	sodium ion transport	18	0	0.08	2.0e-06
GO:0006006	glucose metabolic process	53	0	0.25	6.4e-06
GO:0044267	cellular protein metabolic process	35	0	0.16	2.5e-05
GO:0042026	protein refolding	19	0	0.09	0.00067
GO:0005980	glycogen catabolic process	4	0	0.02	0.00150
GO:0006555	methionine metabolic process	4	0	0.02	0.00462
GO:0006177	GMP biosynthetic process	3	0	0.01	0.00569

for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac.enr, term, ranks=TRUE))
}

## [1] "cellular amino acid metabolic process"

## [1] "oxidation-reduction process"

## [1] "proteolysis"

## [1] "sodium ion transport"

## [1] "glucose metabolic process"

## [1] "cellular protein metabolic process"

## [1] "protein refolding"

## [1] "glycogen catabolic process"

## [1] "methionine metabolic process"

## [1] "GMP biosynthetic process"

#showGroupDensity(GOdata.bac.enr, "GO:0006508", ranks = TRUE) #proteolysis

RAG vs TCell - BioGlyCMK

ontology <- 'MF'

MF

geneList.bac.enr <- -1*geneList.bac.enr 

GOdata.bac.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.bac.enr, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.enr.ks <- runTest(GOdata.bac.enr, statistic = 'ks')
a <- GenTable(GOdata.bac.enr, result.bac.enr.ks, topNodes = 10)

kable(a)

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0050242	pyruvate, phosphate dikinase activity	55	0	0.65	1.3e-25
GO:0016301	kinase activity	61	0	0.72	3.2e-19
GO:0016903	oxidoreductase activity, acting on the a…	141	0	1.66	9.1e-13
GO:0051536	iron-sulfur cluster binding	45	0	0.53	2.7e-11
GO:0016887	ATPase activity	33	0	0.39	9.0e-11
GO:0005524	ATP binding	149	0	1.75	4.7e-10
GO:0004634	phosphopyruvate hydratase activity	44	0	0.52	2.2e-08
GO:0016820	hydrolase activity, acting on acid anhyd…	16	0	0.19	6.8e-07
GO:0000287	magnesium ion binding	57	1	0.67	1.0e-05
GO:0016747	transferase activity, transferring acyl …	26	0	0.31	1.6e-05

for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac.enr, term, ranks=TRUE))
}

## [1] "pyruvate, phosphate dikinase activity"

## [1] "kinase activity"

## [1] "oxidoreductase activity, acting on the a..."

## [1] "iron-sulfur cluster binding"

## [1] "ATPase activity"

## [1] "ATP binding"

## [1] "phosphopyruvate hydratase activity"

## [1] "hydrolase activity, acting on acid anhyd..."

## [1] "magnesium ion binding"

## [1] "transferase activity, transferring acyl ..."

BP - RAG vs Tcell: BioGlyCMK

ontology <- 'BP'
GOdata.bac.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.bac.enr, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.enr.ks <- runTest(GOdata.bac.enr, statistic = 'ks')
a <- GenTable(GOdata.bac.enr, result.bac.enr.ks, topNodes = 10)

kable(a)

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0006090	pyruvate metabolic process	132	4	1.54	1.2e-25
GO:0016310	phosphorylation	133	4	1.55	4.4e-25
GO:0022900	electron transport chain	16	0	0.19	2.9e-05
GO:0005978	glycogen biosynthetic process	4	0	0.05	0.00044
GO:0006096	glycolytic process	77	4	0.90	0.00046
GO:0016226	iron-sulfur cluster assembly	12	0	0.14	0.00054
GO:0006730	one-carbon metabolic process	7	0	0.08	0.00286
GO:0006810	transport	51	0	0.60	0.00352
GO:0005975	carbohydrate metabolic process	224	4	2.61	0.00414
GO:0006091	generation of precursor metabolites and …	116	4	1.35	0.00616

for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac.enr, term, ranks=TRUE))
}

## [1] "pyruvate metabolic process"

## [1] "phosphorylation"

## [1] "electron transport chain"

## [1] "glycogen biosynthetic process"

## [1] "glycolytic process"

## [1] "iron-sulfur cluster assembly"

## [1] "one-carbon metabolic process"

## [1] "transport"

## [1] "carbohydrate metabolic process"

## [1] "generation of precursor metabolites and ..."

#showGroupDensity(GOdata.bac.enr, "GO:0006091", ranks=TRUE) #generation of precursor metabolites
#showGroupDensity(GOdata.bac.enr, "GO:0006090", ranks=TRUE) # pyruvate metabolic process

Enriched vs Unenriched

groups <- read.csv(file.path(data.dir, 'groups.csv'))
groups.rt <- groups[groups$RT_Enriched == 'True' | groups$RT_Unenriched == 'True', ]
geneList.rt.enr <- groups.rt$RT_Enriched == 'True' & groups.rt$RT_Unenriched == 'False'
geneList.rt.enr <- as.numeric(geneList.rt.enr)
names(geneList.rt.enr) <- groups.rt$X

sigFunc <- function(score) {return(score > 0.5)}

MF - Enr vs Unenr: Tcell

ontology <- 'MF'

GOdata.rt.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.rt.enr, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rt.enr.fisher <- runTest(GOdata.rt.enr, statistic = 'fisher')

kable(GenTable(GOdata.rt.enr, result.rt.enr.fisher, topNodes = 10))

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0016620	oxidoreductase activity, acting on the a…	146	76	30.70	1.2e-14
GO:0004816	asparagine-tRNA ligase activity	18	17	3.78	3.8e-11
GO:0003824	catalytic activity	2849	644	599.07	5.0e-08
GO:0008234	cysteine-type peptidase activity	23	19	4.84	3.4e-06
GO:0046872	metal ion binding	472	83	99.25	1.3e-05
GO:0008774	acetaldehyde dehydrogenase (acetylating)…	17	12	3.57	1.5e-05
GO:0004197	cysteine-type endopeptidase activity	15	11	3.15	2.0e-05
GO:0004022	alcohol dehydrogenase (NAD) activity	18	12	3.78	3.6e-05
GO:0004375	glycine dehydrogenase (decarboxylating) …	6	6	1.26	8.5e-05
GO:0004825	methionine-tRNA ligase activity	12	9	2.52	9.1e-05

BP - Enr vs Unenr: Tcell

ontology <- 'BP'

GOdata.rt.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.rt.enr, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rt.enr.fisher <- runTest(GOdata.rt.enr, statistic = 'fisher')

kable(GenTable(GOdata.rt.enr, result.rt.enr.fisher, topNodes = 10))

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0006421	asparaginyl-tRNA aminoacylation	18	17	3.71	2.7e-11
GO:0008152	metabolic process	2834	603	583.98	1.0e-10
GO:0015976	carbon utilization	17	12	3.50	1.2e-05
GO:0006508	proteolysis	89	36	18.34	3.9e-05
GO:0006431	methionyl-tRNA aminoacylation	12	9	2.47	7.7e-05
GO:0006066	alcohol metabolic process	42	13	8.65	0.00014
GO:0009082	branched-chain amino acid biosynthetic p…	24	13	4.95	0.00028
GO:0009081	branched-chain amino acid metabolic proc…	29	18	5.98	0.00034
GO:0006099	tricarboxylic acid cycle	5	5	1.03	0.00037
GO:0006546	glycine catabolic process	8	7	1.65	0.00043

RAG Enr vs Unenr

groups.rag <- groups[groups$RAG_Enriched == 'True' | groups$RAG_Unenriched == 'True', ]
geneList.rag.enr <- groups.rag$RAG_Enriched == 'True' & groups.rag$RAG_Unenriched == 'False'
geneList.rag.enr <- as.numeric(geneList.rag.enr)
names(geneList.rag.enr) <- groups.rag$X

MF - Enr vs Unenr: RAG

ontology <- 'MF'

GOdata.rag.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.rag.enr, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rag.enr.fisher <- runTest(GOdata.rag.enr, statistic = 'fisher')

kable(GenTable(GOdata.rag.enr, result.rag.enr.fisher, topNodes = 10))

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0016620	oxidoreductase activity, acting on the a…	142	77	29.01	4.7e-17
GO:0004816	asparagine-tRNA ligase activity	18	17	3.68	2.3e-11
GO:0003824	catalytic activity	2697	595	550.94	1.2e-08
GO:0008774	acetaldehyde dehydrogenase (acetylating)…	17	12	3.47	1.1e-05
GO:0008234	cysteine-type peptidase activity	22	18	4.49	1.3e-05
GO:0004197	cysteine-type endopeptidase activity	15	11	3.06	1.5e-05
GO:0004022	alcohol dehydrogenase (NAD) activity	18	12	3.68	2.6e-05
GO:0004375	glycine dehydrogenase (decarboxylating) …	6	6	1.23	7.1e-05
GO:0004825	methionine-tRNA ligase activity	12	9	2.45	7.2e-05
GO:0004019	adenylosuccinate synthase activity	8	7	1.63	9.5e-05

BP - Enr vs Unenr: RAG

ontology <- 'BP'

GOdata.rag.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.rag.enr, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rag.enr.fisher <- runTest(GOdata.rag.enr, statistic = 'fisher')

kable(GenTable(GOdata.rag.enr, result.rag.enr.fisher, topNodes = 10))

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0006421	asparaginyl-tRNA aminoacylation	18	17	3.55	1.2e-11
GO:0008152	metabolic process	2674	547	527.01	3.8e-10
GO:0015976	carbon utilization	17	12	3.35	7.3e-06
GO:0006508	proteolysis	84	34	16.56	2.7e-05
GO:0006431	methionyl-tRNA aminoacylation	12	9	2.37	5.3e-05
GO:0055114	oxidation-reduction process	867	207	170.88	6.4e-05
GO:0006066	alcohol metabolic process	42	13	8.28	8.7e-05
GO:0009082	branched-chain amino acid biosynthetic p…	24	13	4.73	0.00017
GO:0009081	branched-chain amino acid metabolic proc…	28	17	5.52	0.00141
GO:0006546	glycine catabolic process	7	6	1.38	0.00146

Unenriched vs Enriched

geneList.rt.un <- groups.rt$RT_Enriched == 'False' & groups.rt$RT_Unenriched == 'True'
geneList.rt.un <- as.numeric(geneList.rt.un)
names(geneList.rt.un) <- groups.rt$X

MF - Unenr vs Enr: Tcell

ontology <- 'MF'

GOdata.rt.un <- new('topGOdata', ontology = ontology, allGenes = geneList.rt.un, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rt.un.fisher <- runTest(GOdata.rt.un, statistic = 'fisher')

kable(GenTable(GOdata.rt.un, result.rt.un.fisher, topNodes = 10))

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0003735	structural constituent of ribosome	153	132	47.98	< 1e-30
GO:0005198	structural molecule activity	298	246	93.44	< 1e-30
GO:0003677	DNA binding	137	111	42.96	< 1e-30
GO:0003723	RNA binding	216	145	67.73	< 1e-30
GO:0032549	ribonucleoside binding	861	279	269.98	< 1e-30
GO:0003899	DNA-directed RNA polymerase activity	137	107	42.96	1.7e-30
GO:0004619	phosphoglycerate mutase activity	44	40	13.80	1.4e-16
GO:0046933	proton-transporting ATP synthase activit…	32	30	10.03	1.4e-13
GO:0016740	transferase activity	869	284	272.49	0.00019
GO:0003924	GTPase activity	119	55	37.31	0.00039

BP - Unenr vs Enr: Tcell

ontology <- 'BP'

GOdata.rt.un <- new('topGOdata', ontology = ontology, allGenes = geneList.rt.un, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rt.un.fisher <- runTest(GOdata.rt.un, statistic = 'fisher')

kable(GenTable(GOdata.rt.un, result.rt.un.fisher, topNodes = 10))

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0071973	bacterial-type flagellum-dependent cell …	112	112	36.01	< 1e-30
GO:0006412	translation	261	173	83.92	< 1e-30
GO:0006351	transcription, DNA-templated	152	116	48.87	9.1e-30
GO:0015986	ATP synthesis coupled proton transport	34	30	10.93	1.3e-11
GO:0000162	tryptophan biosynthetic process	5	5	1.61	0.0034
GO:0006438	valyl-tRNA aminoacylation	5	5	1.61	0.0034
GO:0006556	S-adenosylmethionine biosynthetic proces…	84	39	27.01	0.0039
GO:0015991	ATP hydrolysis coupled proton transport	37	20	11.90	0.0045
GO:0006801	superoxide metabolic process	7	6	2.25	0.0056
GO:0006563	L-serine metabolic process	15	11	4.82	0.0066

RAG Unenriched vs Enriched

geneList.rag.un <- groups.rag$RAG_Enriched == 'False' & groups.rag$RAG_Unenriched == 'True'
geneList.rag.un <- as.numeric(geneList.rag.un)
names(geneList.rag.un) <- groups.rag$X

MF - Unenr vs Enr: RAG

ontology <- 'MF'

GOdata.rag.un <- new('topGOdata', ontology = ontology, allGenes = geneList.rag.un, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rag.un.fisher <- runTest(GOdata.rag.un, statistic = 'fisher')

kable(GenTable(GOdata.rag.un, result.rag.un.fisher, topNodes = 10))

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0005198	structural molecule activity	271	244	92.18	< 1e-30
GO:0003735	structural constituent of ribosome	139	124	47.28	< 1e-30
GO:0003677	DNA binding	136	111	46.26	< 1e-30
GO:0032549	ribonucleoside binding	835	285	284.04	2.2e-30
GO:0003723	RNA binding	202	146	68.71	8.6e-29
GO:0003899	DNA-directed RNA polymerase activity	136	107	46.26	1.3e-27
GO:0004619	phosphoglycerate mutase activity	43	42	14.63	3.7e-19
GO:0046933	proton-transporting ATP synthase activit…	33	31	11.23	5.6e-13
GO:0016661	oxidoreductase activity, acting on other…	23	22	7.82	6.8e-10
GO:0003924	GTPase activity	113	62	38.44	3.0e-06

BP - Unenr vs Enr: RAG

ontology <- 'BP'

GOdata.rag.un <- new('topGOdata', ontology = ontology, allGenes = geneList.rag.un, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rag.un.fisher <- runTest(GOdata.rag.un, statistic = 'fisher')

kable(GenTable(GOdata.rag.un, result.rag.un.fisher, topNodes = 10))

GO.ID	Term	Annotated	Significant	Expected	result1
GO:0071973	bacterial-type flagellum-dependent cell …	112	112	38.53	< 1e-30
GO:0006412	translation	237	167	81.54	< 1e-30
GO:0006351	transcription, DNA-templated	149	114	51.26	2.2e-27
GO:0015986	ATP synthesis coupled proton transport	35	31	12.04	3.4e-11
GO:0015991	ATP hydrolysis coupled proton transport	32	21	11.01	0.00029
GO:0006414	translational elongation	56	31	19.27	0.00094
GO:0000162	tryptophan biosynthetic process	5	5	1.72	0.00479
GO:0006801	superoxide metabolic process	7	6	2.41	0.00812
GO:0006563	L-serine metabolic process	15	11	5.16	0.01086
GO:0006556	S-adenosylmethionine biosynthetic proces…	84	39	28.90	0.01385